Poniżej przedstawiam kod MATLABa obliczający problem adwekcji-dyfuzji metodą Galerkina, dla modelowego problemu Erikksona-Johnsona.
Przypomnijmy najpierw problem adwekcji-dyfuzji, problem modelowy Erikksona-Johnsona (opisany na przykład w dokumencie [1]) zdefiniowany na obszarze kwadratowym \( \Omega = [0,1]^2 \) w następujący sposób: Szukamy funkcji \( u \) takiej że
\( a(u,v)=l(v) \forall v \) gdzie
\( a(u,v) =\int_{\Omega} \beta_x(x,y) \frac{\partial u(x,y) }{\partial x } dxdy + \int_{\Omega} \beta_y(x,y) \frac{\partial u(x,y)}{\partial y } dxdy \\ +\epsilon \int_{\Omega} \frac{\partial u(x,y) }{\partial x} \frac{\partial v(x,y)}{\partial x } dxdy +\epsilon \int_{\Omega} \frac{\partial u(x,y)}{\partial y } \frac{\partial v(x,y)}{\partial y } dxdy \)
\( l(v) = \int_{\partial \Omega } f(x,y) v dxdy \)
\( f(x,y)=sin(\pi y)(1-x) \) jest rozszerzeniem warunku brzegowego Dirichleta na cały obszar, natomiast
\( \beta = (1,0) \) reprezentuje prędkość wiatru wiejącego z lewej strony na prawą, natomiast \( \epsilon = 10^{-2} \) oznacza współczynnik dyfuzji. Wielkość \( Pe=1/ \epsilon = 100 \) nazywa sie liczbą Pecleta, i definiuje ona wrażliwość problemu adwekcji-dyfuzji. Przedstawiona tutaj definicja liczby Pecleta obowiązuje w przypadku jednostkowej prędkości wiatru \( \| \beta \|_2=\sqrt{\beta_x^2+\beta^2_y}=1 \). W ogólnym przypadku rozważa się liczbę Pecleta zdefiniowaną na poszczególnych elementach siatki zgodnie ze wzorem
\( Pe=\| \beta \|_2*h/(2*\epsilon) \). Tak zdefiniowaną liczbę Pecleta użyć można do badania stabilności symulacji, liczba Pecleta mniejsza od jeden daje nam warunek stabilności symulacji na elemencie.

Kod można uruchomić w darmowym środowisku Octave.

Pobierz kod lub zob. Załącznik 3A.

W linii 835 tworzony jest knot vector zbudowany z 5 elementów dla wielomianów kwadratowych
\( knot = simple\_knot(5, 2); \)
W linii 836 ustawiany jest współczynnik dyfuzji (dla współczynnika adwekcji równego 1).
\( epsilon = 2e-2; \)
Kod może zostać uruchomiony w darmowym środowisku Octave.
Kod uruchamia się otwierając go w Octave oraz wpisując komendę
\( advection' \)
Po chwili obliczeń kod otwiera dodatkowe okienko i rysuje w nim rozwiązanie numeryczne oraz dokładne.
Kod oblicza również normę L2 i H1 z rozwiązania
\( Error: L2 24.00 procent, H1 72.08 procent \)
i porównuje do normy L2 i H1 z projekcji rozwiazania dokładnego na przestrzeń funkcji B-spline rozpiętych na danej siatce.
\( Best possible: L2 1.56 procent, H1 53.16 procent \)
Zauważmy że norma te nie będzie równa zero, ponieważ rozwiązanie dokładne nie jest dane funkcjami B-spline, i nawet jeśli znamy rozwiązanie dokładne, to nie jest ona zadane funkcjami B-spline, i jego projekcja na przestrzeń funkcji B-spline również jest przybliżeniem rozwiązania dokładnego. Porównanie to pokazuje jaki bład popełnia nasza metoda numeryczna, i jakie normy mogła by osiągnąć na aktualnej siatce z funkcjami B-spline.
Widać że metoda Galerkina nie daje dokładnego rozwiązania.

Autorzy kodów w MATLABie: Marcin Łoś i Maciej Woźniak.